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SUMMARY 

A computer program was developed for predicting nonlinear uniaxial material 
responses using viscoplastic constitutive models. Four specific models, i.e,, 
those due to Miller, Walker, Krieg-Swearengen-Rohde, and Robinson, are included 
in the present program. Any other unified model can be easily Implemented into 
the program in the form of subroutines. Analysis features include stress- 
strain cycling, creep response, stress relaxation, thermomechanical fatigue 
loop, or any combination of these responses. In this report, an outline is 
given on the theoretical background of uniaxial constitutive models, analysis 
procedure, and numerical integration methods for solving the nonlinear con- 
stitutive equations. In addition, a discussion on the computer program 
Implementation is also given. Finally, seven numerical examples are Included 
to demonstrate the versatility of the computer program developed. 


INTRODUCTION 

The ability to perform accurate structural analysis and design studies, and 
ultimately durability assessments of gas turbine engine hot section components 
depends, to a large extent, on having good material characterization models. 

The material characterization for aerospace applications is made difficult 
because of the complex thermomechanical load history, including elevated tem- 
peratures, through which the components are cycled. The observed material 
deformation phenomena that result from the complex loads and that must be 
predicted accurately and efficiently with these models include the interaction 
of creep and plasticity, cyclic stress-strain hardening (or softening), rate 
dependence, stress relaxation, and creep recovery upon unloading. Generally, 
these phenomena are modelled with macroscopic (continuum) constitutive 
theories, which include classical plasticity theories, unified theories, 
decoupled theories and rheological theories. Some of these theories have been 
used with varying degrees of success, where each theory has its unique desir- 
able features and modeling capabilities. 

Recently, there has been considerable research effort on the development of 
unified theories (viscoplastic constitutive models) to predict the inelastic 
behavior of metals. Because the viscoplastic models can represent the inter- 
action of creep and plasticity (an important material deformation phenomenon 
in engines) they have been applied to the characterization of hot section com- 
ponent materials (refs. 1 and 2). Classical theory of plasticity, on the other 
hand, cannot model this phenomenon, as discussed in reference 3. In classical 
theories the inelastic deformation is partitioned into independent plastic and 
creep terms, whereas in viscoplasticity theory it is combined into a single 



Inelastic strain term. Also, the other material deformation phenomena, such as 
cyclic hardening or softening, rate dependence and creep, can be better repre- 
sented with the viscoplastic models than with the classical and other models. 
Consequently, a number of viscoplastic models have been proposed for more 
realistic representation of material deformation of turbine engine hot section 
components . 

Viscoplastic models, like classical models and others, are not, however, 
without their shortcomings. Some of the shortcomings are: (a) the models 

have not been fully developed or adequately tested; (b) the determination of 
material parameters In these models Is difficult since It Is based largely on 
trial-and-error fitting of experimental data; (c) there Is only a limited 
experimental data base available for materials used In hot section components; 
and (d) the associated constitutive differential equations have stiff regimes 
that present numerical difficulties In time-dependent analyses. To take full 
advantage of these models, the aforementioned deficiencies must be overcome. 

It Is also noted that the differential equations of the viscoplastic models 
are highly nonlinear. Thus analytical solution of these equations, except In 
a few special cases. Is Impossible and numerical solutions are required. 

In this report, a computerized analysis procedure to predict the nonlinear 
uniaxial stress-strain response of materials subjected to a variety of thermo- 
mechanical loads (simple to complex) using several viscoplastic models was 
developed. This procedure was Incorporated Into a computer program, called 
VPMOOEL, with the following purposes In mind: 

1. To compare the responses predicted by a theory with experimental data. 

2. To help understand the physical phenomena represented In a particular 
viscoplastic model . 

3. To perform parametric studies. 

4. To conduct comparative studies of several viscoplastic models. 

5. To Investigate different numerical Integration schemes. 

The computer program Is a modular structured code such that any viscoplastic 
model or numerical Integration scheme can be added or deleted from the program 
without upsetting Its overall coding organization. In addition, all computa- 
tions for uniaxial responses under monotonic, cyclic, creep or stress relax- 
ation loads, or any combination of these loads, are performed In an Interactive 
mode with free-format Input. A computer menu guides the analyst through 
the program, even to graphics plotting options for quick analysis of the data. 

Included In this report Is a section on the background of four visco- 
plastic models. Each model Is briefly reviewed and the differential equations 
are presented In their uniaxial form. Also, Included In other sections are 
discussions of the analysis procedure and various numerical Integration schemes 
for solving the uniaxial stiff differential equations, computer Implementation 
of these equations, and numerical examples to demonstrate the utilization of 
the computer program VPMODEL. 


BACKGROUND 

Numerous viscoplastic models have appeared In the literature. Although 
they differ In details, most of the models share several common features. For 
example, the Inelastic strain rate Is assumed to be a function of stress, 
loading history, and loading rate. Such dependence Is Introduced through the 
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use of two state variables, namely the equilibrium stress a and drag stress 
k. Moreover, the total strain rate is assumed to be the summation of the 
elastic and Inelastic strain rate terms. With these assumptions, a skeletal 
form of vlscoplastlc models can be written as follows: 



where 

ej = inelastic strain 
a = applied stress 
a = equilibrium stress 
k = drag stress 
( •) = time derivative 

h^j, h|( and are the work-hardening and recovery functions 

respectively; they are generally functions of applied stress, temperature and 
state variables. Equation (1) is the inelastic strain rate or flow law and 
equations (2) and (3) are called evolutionary equations for the state vari- 
ables. In addition to the above relations, we invoke the assumption that the 
stress rate is proportional to the elastic strain rate by 


o = E (e - 


•I 

e 


- BT) 


(4) 


where E is the Young's modulus of the material; c is the total strain; B is 
the thermal expansion coefficient; and T is the temperature. 


By assigning appropriate terms in the evolutionary equations, several 
important deformation phenomena of metals can be simulated by the viscoplastic 
models. For example, by allowing the equilibrium stress a and drag stress k 
to vary with an increase in cumulative inelastic strain in equations (2) and 
(3), cyclic hardening (or softening) of hysteresis loops of the material can be 
modelled. The transient creep or stress relaxation phenomena are represented 
by the rapid growth of the state variables during the initial stage of loading. 
Such growth becomes saturated as the recovery terms in equations (2) and (3), 
the second terms, become dominant. That is, steady state creep is reached when 
the strain work-hardening (or softening) process is in dynamic equilibrium with 
the thermal recovery process. Other metallurgical phenomena, such as the 
effects of annealing and warm-working, can be included by introducing appro- 
priate terms in the evolutionary equations. In addition, temperature depend- 
ence of the material can also be included by allowing the inelastic strain and 
state variables to change with temperature. One unique feature of the vlsco- 
plastlc theory, differing from the classical plasticity theory, is that it does 
not Involve the concept of a yield surface explicitly. 


By proposing various mathematical expressions for the inelastic strain and 
state variables, a number of viscoplastic models have emerged in the litera- 
ture. The more notable models are those due to Hart (ref. 4), Miller (ref, 5), 
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Walker (ref. 6), Kr1eg, Swearengen and Rohde (ref. 7), and Robinson (ref. 8). 
Since the Hart model can be considered as a special case of Walker's, it will 
not be discussed herein. In this section, only the latter four models are 
briefly outlined in their uniaxial differential forms. 


Miller's Model 


The inelastic strain rate for Miller's model was assumed to be a hyperbolic 
sine function of the applied] stress and two state variables. The rate equa- 
tions for inelastic strain and the two state variables for Miller's improved 
model (ref. 9) are^ 


• I 

C 


Be' jsinh 


sgn (<j/E - o) 


a 



a = - H.|Be' [sinh (A.,|a|)]'’ sgn (a) 

\a\ - ~ • k^-^j Ic^l Mind.O, e"/e') 


( 5 ) 

( 6 ) 




- H2C2B6" [s 1 nh(A 2 K^-^]" ( 7 ) 


1 

I where 


»1 

= H.j exp [- H2osgn(c^)] 

(8) 


e' 

= exp {(-Q*/RT^) • In(VT) + 1 ]) if T < T^ 

( 9 ) 


e' 

= exp(-Q*/RT) if T > T^ 

(10) 


e" 

= X" exp . [ln(T^/T) + 1 ]| if T < T^ 

(11) 



e" = X" exp {- Q^e^o/RT} T > T^ 

(12) 


with X" 

= exp(RpX') 

( 13 ) 



'Some of the terms in Miller's latest model which simulate metallurgical 
behavior have not been included herein. with 
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(14) 


Q 

recov = Q*[l + (R - 1 ) X' ] 

Q 


2 

X' = (1 - X) 1f X > 1 


X' = 0 1f X < 1 


(15) 


X 


q/E - tt 

2/3 1.5 

(A ) (k) 

2 


(16) 


In the above equations B, E, n, H] , H 2 , H 3 , A], ^ 2 • ^2» Q*» mate- 

rial constants, which, except for E, are temperature Independent; Tt Is the 
transition temperature of metal. Temperature dependence of the model Is 
Introduced through the use of factors o' and e" (activation energies). The 
model considers both kinematic and Isotropic hardening of the material due to 
the growth laws for o and k. In the original Miller's model (ref. 5), the 
equilibrium stress was assumed to harden linearly with Inelastic strain. This 
circumscribes a tri-llnear curve In modeling cyclic stress-strain response. 

The model was subsequently modified (ref. 9) to accommodate nonlinear hardening 
effects by allowing the constant HI to vary with <x according to equation ( 8 ). 


Walker's Model 

This model was derived from a three 3-element mechanical model for visco- 
plastic materials (ref. 6 ). In Its mathematical form, the Inelastic strain Is 
expressed by a power law In terms of the applied stress and two state vari- 
ables. The equations for the Inelastic strain rate and the rate dependent 
state variables are 


•I 

c 


sgn(q - a) 


(17) 


* *I I "”^5lc^l 

O = (n.| + n^)^ - (a - <xq - n.|G ) {(n^ + n^e ') 


, • 1 1 I ,m-l 

|c I + n^kl } 


(18) 
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In the above equations there are nine material constants, namely m, n, n-| , 

^2* f'S* ^4* ^5* f’S* *17* “0* *^1 nonisothermal conditions, all 

material constants may vary with temperature. Since the model contains two 
state variables, I.e., equilibrium stress and drag stress, both the kinematic 
and Isotropic hardening effects can be simulated. The growth law for equilib- 
rium stress (eq. ( 18 )), Includes both dynamic and static thermal recovery 
terms; whereas only the static recovery was considered for the growth of drag 
stress, (eq. ( 19 )). Walker's model was formulated In both differential and 
Integral forms (ref. 1 ). Only the differential form Is considered In this 
report. 


Krieg, Swearengen and Rohde's (KSR) Model 

The KSR model (ref. 7 ) Is similar to Walker's model, since It also uses a 
power law to represent the Inelastic strain rate. The Inelastic strain rate 
equation Is 


•I lo-Oll*^ / V 

c = ^ sgn(o - a) 


while the state variable rate equations are 


a = - ^2<x la| y® ^ ^ 

k = 1^1 - Ag(k - kg)" 


The seven material constants In these equations are n, A^ , A2, A3, A4, 

A5, and kQ, which can be temperature dependent. Both the equilibrium 
stress and drag stress are Included In the formulation (I.e., eq. ( 20 )). The 
growth laws contain two terms, the linear hardening and static recovery terms. 
These terms produce tri-llnear cyclic stress-strain curves or hysteresis loops. 
This behavior may be modified by assuming the constant A3 In equation ( 21 ) 
to be very small. Furthermore, cyclic hardening could be modelled by letting 
the constants A4 and A5 have small values In equation (22) so that the 
drag stress grows slowly with deformation. 


( 20 ) 

( 21 ) 

( 22 ) 


Robinson's Model 

This model (ref. 8), to some extent. Is also analogous to the previous 
ones. That Is, the Inelastic strain rate Is assumed to be a power function of 
the effective stress and two state variables. The equations for the Inelastic 
strain rate and the rate dependent state variables are 
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where 


•I , a-P(2^j^) < _ i>n sgn(a - a) 

a = Hla'T® - R|a"r“® sgn(a) 
k = 0 

“■ ■ - “0> (w^) * “0 


r 


a 3 


a . I<»l > 2aj 

2 

fc: * ‘"0* "" ^“c 


n 


P(x) = { 


1 - (1 - xr/2 
(1 + x)^/2 
LO 


X > 1 
0 < X < 1 
-1 < X < 0 
X < -1 


<x> = 



X > 0 
X < 0 


(23) 

(24) 


(25) 


(26) 


(27) 


(28) 


The eight material constants In these equations are A, n, k, B, m, H, ag 
and R; W-| and W2 are smoothing constants to account for discontinuities In 
the material characterization and their values are selected by the analyst. 

For nonisothermal conditions, temperature dependence of the model Is allowed 
by varying the constants A and R with temperature. In the above equations, 
the drag stress was assumed to be constant. Thus, the current version of 
Robinson's model only considers the kinematic hardening of materials for which 
Its Isotropic hardening has already become saturated. If Isotropic hardening 
Is desired. It can be easily Introduced with a growth law for k In the form 
similar to equation (22) (KSR's model). 


ANALYSIS PROCEDURE 

With the uniaxial viscoplastic models outlined In the previous section. 
It Is possible to analyze three basic types of material response. The types 
of response are 

(1) Monotonic or cyclic stress-strain response 

(2) Creep response 

(3) Stress relaxation 
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Also, any combination of the above responses can be obtained by specifying 
appropriate loading histories. The analysis may be performed for either 
Isothermal or nonisothermal conditions. Caution should be given, however, to 
nonisothermal analysis, since the constitutive model must be properly defined 
for such a condition. Each of the basic response types are briefly described 
below. 


Monotonic or Cyclic Stress-Strain Response 

For viscoplastic materials, time dependency Is an Important consideration 
In the calculation of material response. As an example, for the purpose of 
describing the analysis procedure, consider for a given strain rate the cyclic 
stress-strain response of a material. In this case It Is divided Into several 
segments for the purpose of proceeding with the analysis: (a) loading; (b) 

unloading (or reverse loading); (c) loading again, etc. For each loading (or 
unloading) segment, the analyst must supply the following Information to begin 
the analysis: 

Temperature T 

Strain rate c 

Time Increment at 

Strain limit ci 

Using the strain limit, a stress limit oi may be supplied to terminate the 
analysis of a load segment. Based on this Information, the following steps 
are used In the analysis procedure to calculate the stress and strain at any 
time tn+i : 


Step 1. Calculate the Increments of total strain and Inelastic 
strain (assuming Isothermal condition for discussion 
purposes) according to 


Ac s c • At 

A I A^- 

Ac = e • At 


(29) 


where c^ Is computed from equation (1). The specific expression of c^ 1s 
obtained from whichever viscoplastic model Is selected. 


Step 2. Calculate the stress rate and stress Increment from 


• • • I 

= E(c - c ) 

• 

Acf s a • At 


(30) 
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step 3. Update the total stress, total strain and the state 
variables Involved. 


c(t ,) = c(t ) + Ac 

n+r ' n' 


a( t - ) = a( t ) + Ao 
' n+1 ' n' 


( 31 ) 


= k(t„) * «k 


The above calculation steps are repeated until the total strain or stress 
reaches the specified limit. 


Creep Response 

Creep response Is generally defined as the time-dependent deformation of 
a material under constant stress. In the laboratory, for example, a uniaxial 
stress Is applied to a specimen In a two-step loading history to simulate 
creep. That Is, the stress Is Increased at a constant rate to a prescribed 
value In a short time Interval, say [0, t-|], then held constant afterwards 
as shown In figure 1(a), where ti « t 2 * In addition to two-step load- 
ing, creep tests may also be conducted under multiple-step loading as shown In 
figure 1(b) or under constant stress at any point on the loading/unloading 
branch of hysteresis loops (fig. 1(c)). To consider these different cases, a 
creep analysis can be conveniently performed by dividing Its loading history 
Into two or more segments. For example, the two-step loading shown In figure 
1(a) Is divided Into: segment 1 - stress Is Increased form zero to a constant 
value, and segment 2 - stress Is held constant. The required data the analyst 
must supply for each of the two load segments are 

Segment 1 

(a) Temperature T 

(b) Strain Rate c 

(c) Time step size At 

(d) Strain limit cl or stress limit <jl 

Segment 2 

(a) Temperature T 

(b) Time step size At 

(c) Hold time t^ 

Division of the loading history for multistep or cyclic loadings can be made 
by repeating this procedure. 


9 



Stress Relaxation 


In contrast to creep, stress relaxation Is obtained by applying a con- 
stant strain to a specimen and holding this strain for an extended time 
period. In this case, the analyst must specify: temperature, strain rate, 

time Increment, strain limit and hold time. The analysis procedure Involved 
herein for two-step or multlstep loadings Is almost Identical to that of creep 
analysis. 


NUMERICAL INTEGRATION 

In order to calculate the uniaxial response (either stress or strain) of 
a material using a viscoplastic theory, a system of highly nonlinear differen- 
tial equations, such as those of equations (1) to (4) must be solved, by 
numerical Integration. It Is well known that these equations have "stiff" 
regimes and special attention Is, therefore, required to solve them numeri- 
cally. For example, the Inelastic strain rate In equation (1) Is a strong 
nonlinear function of the stress and the state variables a and k. That Is 
to say, any small variation In o, a and k can cause significant changes In 

the value of c^. Similar behavior Is also found In the growth laws (I.e., 
eqs. (2) and (3)). Therefore, a very small time step Is often required to 
Integrate these equations by numerical methods. 

In this section, various methods for solving the nonlinear differential 
equations of viscoplastic models are discussed. The question of convergence 
control and possible ways of selecting the size of the time step are also con- 
sidered. For the purpose of discussion, equations (1) to (3) are replaced by 
the following matrix expression 


i = t(x.t) 


(32) 


This expression represents a system of nonlinear ordinary differential equa- 
tions. Although various numerical schemes can be used to solve these equa- 
tions (e.g., refs. 10 and 11), practical considerations must be given to three 
Important Issues: (a) their suitability for large scale stress analysis; (b) 

their solution accuracy; and (c) computer solution time. For Instance, Gear's 
multlstep method (ref. 12) was utilized effectively by Miller (ref. 13) for 
obtaining uniaxial viscoplastic response of Zircaloy, but the method Is not 
suitable for Incorporation Into a finite element analysis due to the coupling 
of global degrees of freedom for multiple time steps. Moreover, the method 
requires an Independent procedure to start-up an analysis (ref. 14). 

In view of the above discussion, three simple numerical schemes were 
selected for the present analysis, namely the forward Euler (explicit), 
explicit trapezoidal and Implicit trapezoidal Integration schemes. 

In the forward Euler Integration scheme, the value of y at time tn+i 
In equation (32) Is approximated by 
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( 33 ) 


t 


where 


^n+1 “ * ~i 

^n+1 * 

in ■ J!‘tt,> 


If the above equation 1s employed In an analysis of the "stiff* differential 
equations, very small time steps must be enforced to avoid any numerical 
Instability. 


Instead of using equation (33), the value of y at time tn^i may be 
approximated by an Implicit trapezoidal rule 


^n+1 " ^n ^ 2 


At ,, 

^ O (f,f 


~n+l ^ 


(34) 


where f^+i = f (yn+i)* With the application of the Newton-Raphson method, 
equation (34) Is rewritten for the 1-th Iteration 


[I - At . Ay^^^ = F 


Un+i Xn+1 


at) 


(35) 


where I = Identity matrix 

af 


J = — (A nonsymmetric Jacobian matric) 

Oj 


- Xn 


(36) 


F = V + at«f^/2 


(37) 


(Eq. (35) Is also called an Implicit trapezoidal scheme with Iterations.) If 
this method Is employed In an analysis, the Immediate question Is - how can It 
be determined whether the solution has converged or not? Several convergence 


criteria could be used for this purpose. 
Iterative value of Ay such that 

e . 

II I II 


One convenient way Is to check the 
< Tol. (38) 


where 


II II = Euclidean norm 


Tol = Tolerance ratio 

Presently, the above criterion Is Implemented Into the computer program VPMODEL 
to determine the convergence of a solution. The Implicitness of f^+i In 
equation (34) may be removed by approximating 
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~n+l 

= f + 
~n 

Vi 

(39) 

•’n 

\3y/ 

C 

II 

(40) 

Thus, equation (34) becomes 




[i- 

^ J 1 

2 «n 

Ay = At • f^ 

(41) 


The above equation 1s called the explicit trapezoidal scheme. These three 
schemes have been Implemented Into VPMOOEL for the study of numerical efficiency. 

Another concern Is the selection of the time step. It Is possible to 
estimate a time Increment If the strain rate Is expressed by a power law 
(Norton's law) of the form 


e = A o" (42) 

where A and n are material constants and 

Og = (o - a)/k (43) 

When the Euler forward method Is employed, a step size for obtaining a stable 
solution Is given by (ref. 15) 


At 


s 



(44) 


where 


Atj = Largest time step that can be used for obtaining a stable 
solution 


and 


E = Young's modulus 


Obviously, the above estimate Is valid only If the elastic strain rate Is much 
smaller than the Inelastic rate. If the constitutive equation for Inelastic 
strain rate Is considerably different from the Norton's power law, then an 
approximate relationship In the form of equation (42) can be established. 



with 1 _ d Ing 

n ■ d In* 


(45) 

(46) 
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In view of the above discussion and computer experimentation runs using 
VPMOOEL, the following step size Is suggested for high strain-rate 
(~10~2 _ iO~^/sec) monotonic or cyclic loading: 

Euler forward Integration scheme At = atj 

Trapezoidal Integration with or without Iterations At = SAtj 

For steady-state creep or stress relaxation a much larger step size may be 
employed. 

For different types of loading, an automatic time step control Is useful, 
first to speed-up computation time and second to reduce data storage require- 
ments. This Is done according to the following criteria: 

Let 

At a trial step size 
e calculated error 
®max maximum error limit 
®m1n minimum error limit 
At step size used for analysis 

then 

(a) e > emax» t = At/2 

(At Is further halved until e < e^gx) 

(b) ejnin < e £ ©max* = At 


(c) e < einin» ^t = 2At 

(This process Is repeated untll^emin < ® £ ®max*) 

The definition of the error e varies depending on which Integration scheme Is 
to be used. The error Is measured by 

e ~ vff, for the Euler scheme 

and 

e ~ for the trapezoidal scheme 


where 


= ^n “ ^n-1 


= fn - 2fn-l + ^ 


(47) 

(48) 


These equations hold for the backward difference approximation. In view of 
the above relationships, the measures for error are chosen to be 
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2 

11 7 fn ii *At for the Euler scheme 
® “ II yn II 

2 

liTfnii »At for the trapezoidal scheme 
® “ II yn II 


(49) 

(50) 


On the other hand, 1f the Implicit trapezoidal scheme Is employed, It Is 
more convenient to adopt the relationship given In equation (38) for the 
calculation of error. 


COMPUTER IMPLEMENTATION 

The analysis procedure described has been programmed and the four visco- 
plastic models discussed have been Incorporated Into a computer program called 
VPMOOEL to calculate the uniaxial responses of materials. The uniaxial 
responses Include cyclic stress-strain history, creep, stress relaxation, and 
any combination of these phenomena. It Is Important to recall that the pur- 
poses of writing this program are Intended for: (a) correlating the responses 

predicted by a theory with experimentally measured data; (b) performing para- 
metric studies of material constants for a particular viscoplastic model, (c) 
performing comparative studies of various viscoplastic models; and (d) Investi- 
gating different numerical Integration schemes In terms of their computational 
efficiency. 

With the above purposes In mind, the following features are Included In 
the program: 

(1.) The viscoplastic models and Integration schemes were coded (In 
Fortran) In separate modules so that they can be easily modified, added to or 
deleted from the program. 

(2.) All computations are In an Interactive mode, and the Input data are 
specified In conversational form and free-fleld format. 

(3.) Analysis results can be plotted graphically and/or printed 
numerically. 

At present, material constants of two high temperature alloys for the specific 
models listed below are Included In the program: 

Hastelloy X - material - Miller's, Walker's and KSR's models 

2-1/4 Cr-1 Mo Steel - Robinson's model 

Other material types or material/model combinations. If known, may be added to 
the program quite easily by the user. 

The program begins with four primary decisions to be made by the user, 
namely 

(1) Material model selection 

(2) Material type 
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(3) Integration scheme to be used 

(4) Analysis type 

A menu for each of the above four items is listed in figure 2. After the 
user has made the selections, VPMOOEL will proceed with the corresponding 
analysis in accordance with the analysis procedure outlined In section 3. At 
the end of each load segment, the user has several options to either continue 
the analysis, plot the results, or terminate the execution. A flow chart for 
the analysis of VPMOOEL is shown in figure 2. 

The graphics part of the program was written by using a Calcomp plotting 
package with the plotting device being either the Calcomp plotter or the 
Tektronix graphics terminal, or their equivalent. The free-field input format 
was coded externally so that the program does not rely on the availability of 
such software on a particular compiler. 


EXAMPLES 

To demonstrate the utility of the program VPMOOEL, seven example problems 
are included in this section. Two alloys are considered, namely Hastelloy-X 
and 2-1/4 Cr-1 Mo ferretic steel. Hastelloy-X is used for jet engine combustor 
liners whereas the ferretic steel is a typical steam generator material for ad- 
vanced nuclear reactors. The material constants of Hastelloy-X for the Miller, 
KSR and Walker models are given in references 1 and 2. Material constants for 
the ferretic steel with Robinson's model are reported in reference 8. 


Cyclic Stress-Strain Response Using Walker's Model 

The cyclic response of Hastelloy-X, represented by Walker's model, is 
considered first. In this case, the following analysis parameters were 
specified: 

Temperature T = 1600® F (871® C) 

Strain rate c = 3.87 x 10-3/sec 

Strain limit gl = 0.006 

Step size At = 0.04 sec 

Integration method = Explicit trapezoidal scheme 

At first, one full cycle of load, consisting of three loading segments, was 
imposed as shown below 

Segment 1 - Loading 0 < e < gl 

Segment 2 - Unloading cl ^ i - ®L 

Segment 3 - Loading -gl < c < gl 

Corresponding to the above three loading segments, the following response 
curves were plotted: 

(a) Stress versus Strain (fig. 3) 
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(b) Back stress versus Inelastic Strain (fig. 4) 

(c) Back stress versus Total Strain (fig. 5) 

To obtain the cyclic stress-strain response of the material, another unloading 
segment (segment 4; -cl < c < cl) was added to the previous three segments 
and the response was plotted for the last two loading segments as shown in 
figure 6. The purpose of applying two extra loading segments is to ensure 
that Hastelloy-X has reached its saturated state due to cyclic hardening. 
Similar responses were obtained for different strain rates at the same 
constant temperature, i.e.. 

Strain rates of = 3.66x10“^, 3,70x10“^, 1.11x10“^ and 1.25x10"^ 

1/sec and corresponding step sizes of t = 0.4 sec, 4 sec., 12 sec, 80 
sec. 

The above results are also shown in figure 6. These curves are almost iden- 
tical to those in figure 59 of reference 1. 


Stress Relaxation Test Using Walker's Model 

A stress relaxation case can be obtained by imposing a constant strain 
history to the material concerned. The relaxation phenomenon of Hastelloy-X 
based on Walker's model is illustrated with the following loading condition: 


Constant temperature 

Strain rate 

Strain limit 
Hold time 
Initial step size 


T = 1600“ F (871“ C) 

c = 3.87 X 10-3/sec 

CL = 0.006 

tn = 50 sec. 
at =0.40 sec. 


Initially, the material was subjected to one full cycle of strains, -cl < c 
< CL (see fig. 7). Then it experienced stress relaxation under constant 
strain (c = cl) with a hold time tn = 50 sec at the maximum stress point 
as shown in figure 7. The total load history was applied in four segments: 


Segment 1 
Segment 2 
Segment 3 
Segment 4 


Loading 0 £ e £ «l 

Unloading cL ^ ® - “®L 

Reloading -cl < c £ cl 

Constant Strain c = cl 


The integration scheme used was the implicit trapezoid rule with variable step 
size for all the four load segments. The calculated stress history (or relax- 
ation curve) is shown in figure 8, where the stress history corresponding to 
the first strain cycle is omitted. At the beginning of hold time, the stress 
value is about 37 ksi (255.3 MPa) and it reduces to 12.5 ksi (86.3 MPa) asymp- 
totically after undergoing 50 seconds of relaxation. 
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Cyclic Stress-Strain Response Using Robinson's Model 


* 


f 


The third example 1s to obtain cyclic stress-strain response of 2-1/4 


Cr-l Mo steel using Robinson's model, 
stress with the following conditions: 


Constant temperature T = 

Strain rate e = 

Strain limit cl = 

Step size At = 

Integration scheme 


The material is subjected to uniaxial 


1000“ F (538“ C) 


2.4 X 10-2/hr 

0.0032 
2 X 10-5 hr 

Explicit trapezoidal method 


Initially, the uniaxial loading was applied in 3 segments: (a) loading; 

(b) unloading; and (c) reloading, similar to example 1. Corresponding to the 
above loading and unloading segments, two response curves were plotted. 

(a) Stress versus Strain (fig. 9) 

(b) Back Stress versus Inelastic Strain (fig. 10) 


It is seen in figure 9 that corresponding to the same strain limit, the stress 
for the reloading segment is higher than the stress of the initial segment 
(virgin material). This phenomenon is associated with strain hardening effect 
portrayed by the model . The hardening effect becomes saturated after the 
material has undergone one full cycle of loading. Therefore, in order to 
obtain a stable cyclic stress-strain response, another unloading segment was 
added and the corresponding stress-strain plot is shown in figure 11. In 
addition to the strain rate of 2.4xlQ-2/hr, two different strain rates were 
analyzed, i.e., = 0.24x10-2 and 0.024xl0-2/hr , and the results are plotted 
in figure 11 . 


Creep Response Using Robinson's Model 


To demonstrate how to obtain the creep response by using VPMODEL, we 
consider Robinson's model as an example. The test condition is specified as 
follows: 


Constant temperature T 


1000“ F (538“ C) 


Strain rate 
Stress limits 

Hold times 


c s 0.6xl0-3/hr 

ai = 8, 10, 11, 12.5, 15 ksi 

(ksi = 6.9 MPa) 

tn = 2000, 2000, 2000, 870, 200 hr 


Before the creep loading, a full cycle of stress, o = ± ol» was applied to 
the material so that its strain hardening became saturated. Then the stress 
was held constant at the maximum tensile stress point, i.e., o = ol» 
designated hold time. In the analysis, the implicit trapezoidal scheme with 
variable step sizes was employed. The initial step size for all stress values 
was: At = 1 hr, then the program computed required step size according to 
internally specified convergence tolerance. The calculated creep responses 


17 


(Inelastic strain versus time) corresponding to the hold times are shown in 
figure 12. Again, the creep curves obtained from VPMODEL are identical to 
those generated by Robinson (ref. 8). 


Comparison of Three Viscoplastic Models 

In this example, the stress-strain response of Hastelloy-X material is 
predicted by the three different models (Walker, KSR and Miller) for the same 
temperature, strain rate and loading condition. For this purpose, we have 
considered five cases by varying the temperature and strain rate as defined in 
table 1. Only monotonic loading (one load segment) was imposed for each case 
and the stress-strain plots for all five cases are shown in figures 13 to 17. 
Ideally, for the same material under the same loading condition, the predicted 
responses from the different viscoplastic models should be the same. This is 
not the case according to the results shown in figures 13 to 17. If we use 
Walker's model as a basis for comparison, KSR model correlated with Walker's 

prediction quite well for T = 1400® F (760® C), c = 3.87xl0-3/sec, whereas 

Miller's model under-estimates the stress after the inelastic strain is initi- 
ated as shown in figure 13. On the other hand, as shown in figures 14 and 15, 
for the same temperature but at intermediate or lower strain rates. Miller's 
model correlates more closely with Walker's model than KRS model. Similar 
qualitative results are shown in figures 16 and 17 for a higher temperature, 

T = 1800® F (982® C). 


Numerical Stability Using Walker's Model 

To illustrate the numerical instability problem, a monotonic load case is 
considered. The stress-strain response based on Walker's model was calculated 
for three different step sizes using the forward Euler Integration scheme, 
with t = 0.08, 0.32 and 0.40 sec. The results are shown in figure 18. It is 
obvious that significant numerical oscillation is seen for t * 0.40 sec and 
the oscillation diminishes as the step size is reduced. To make a further 
qualitative comparison of the three integration methods in VPMODEL, we 
consider again the same stress-strain response of Hastelloy-X using Walker's 
model. The problem was run by using the Euler forward, explicit and implicit 
trapezoidal methods (no interactions). According to the results shown in 
figure 19, the trapezoidal methods are far more accurate than the Euler method. 


Thermomechanical Fatigue Loops 

Thermomechanical fatigue loops (TMF) are defined in this study to be 
typical (or simulated) loading histories of turbine engine combustor liners. 
Under these realistic loading conditions, both mechanical load (in the form of 
imposed strain) and temperature are subjected to large changes as functions of 
time. Having accurate predictions of the stress/strain response for these 
loading conditions is essential in the structural integrity and durability 
assessments of combustor liners, particularly at the critical locations in a 
combustor liner (usually where fatigue occurs). To this end, research programs 
have been conducted to develop time-temperature dependent constitutive 
relationships and improve finite element structural analysis procedures 
(refs. 1 and 2). 
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Rather than using sophisticated finite element structural analysis 
programs for evaluation and comparison of the viscoplastic models for dif- 
ferent TMF loop, the uniaxial stress/straIn responses of Hastelloy-x can be 
determined from the viscoplastic models In VPMOOEL. Considered herein for 
Illustrative purposes are three TMF loops: Case 1 - Closed symmetric loop. 

Case 2 - Open symmetric loop, and Case 3 - Open nonsymmetric loop. With these 
specified strain and temperature histories, stress responses were calculated In 
VPMOOEL with the Walker, KSR and Miller models. These examples are considered 
for two reasons: (1) to demonstrate the capability of VPMOOEL In handling non- 

Isothermal loadings; and (2) to compare the predicted uniaxial responses with 
the experimental data reported In reference 2. 

For Case 1, the temperature varies sinusoidally with time from 800* to 
1600* F (427* to 871* C) with a period of 1 min (see fig. 20). Plots of the 
corresponding strain histories (also varying In a sinusoidal form) and the 
strain-temperature relationship are shown In figures 21 and 22. Analysis 
results obtained from the three models In VPMOOEL, which are Identical with 
those from reference 2, are compared with experimental data In figures 23 to 
25. It Is noted that the experimental data stabilized after 2 cycles. 

Figure 23 shows that prediction made by Walker's model requires several cycles 
of loading-unloading before the response stabilizes, and the model considerably 
over-predicts the stress response. The prediction made by Miller's model, 
shown In figure 24, tends to reach the steady-state In the third cycle; 
however, the stress-strain response becomes almost a straight line as It 
approaches the steady-state condition. In other words, the hysteresis effect 
Is erased by the effect of temperature, which Is contrary to the experimental 
data. The prediction made by the KSR model, on the other hand, shows very 
good agreement with the experimental measurements as seen In figure 25. 

For the case of the open symmetric loop. Case 2, the temperature history 
(fig. 26) Is Identical to that of Case 1: The strain history and strain- 

temperature plots are given In figures 27 and 28. Comparisons between 
predictions and experimental data are shown In figures 29 to 31. Both the 
Miller and KSR models attain the steady-state condition during the second TMF 
loop, while Walker's model stabilizes In the fourth cycle. Tensile yields at 
two locations In the loading branch of the hysteresis loops, as Indicated In 
figures 29 and 31, were predicted by both the Walker and KSR models, while 
Miller's model did not reflect such an effect (see fig. 30). For this loading 
condition, the predictions by Walker's and Miller's models are reasonably close 
to the experimental data, while the KSR model over-predicted stress reduction 
for the reverse loading branch of the TMF loop. 

The temperature and strain histories, and strain-temperature relationship 
for Case 3 are shown In figures 32 to 34, respectively. In this case, the 
temperature varies sinusoidally from 950* (510* C) and 1750* F (954* C) for 
1 min with a hold time of 40 sec at the peak temperature. The strain Is held 
constant at -0.43<for the hold-time period. Analysis results for the three 
models and experimental data are compared In figures 35 to 37. It Is seen 
from these figures that In approaching the steady-state condition, the same 
behavior Is exhibited as the previous loading case. Since the constant strain 
hold time corresponds to a constant temperature, theoretical predictions 
should show some degree of stress relaxation. This Is Indeed the case as seen 
in the stress-strain plots (figs. 35 to 37). However, the experimental data 
do not show the distinct stress relaxation changes as do the predictions. 
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Nevertheless, stress relaxation may have occurred gradually during the experi- 
ment. Miller's model predicts a larger amount of stress relaxation than the 
Walker and KSR models. Furthermore, Miller's model does not produce tensile 
yield In the tensile portion of the cycle and simply exhibits elastic behavior 
as shown In figure 36. Because of these two factors, the predictions due to 
the Walker and KSR models are slightly better than those from Miller's model 
for the third loading case. 


CONCLUSIONS 

A computer program called VPMODEL was developed to calculate the uniaxial 
responses of material by using four viscoplastic models, namely Miller, Walker, 
Krieg-Swearengen-Rohde, and Robinson models. The nonlinear uniaxial responses 
which can be analyzed are cyclic stress-strain, creep, stress relaxation, 
thermomechanical fatigue loop or any combination of these phenomena. In 
addition, three Integration schemes, I.e., Euler forward difference, explicit 
and Implicit trapezoidal methods with Newton-Raphson Iterations, were Imple- 
mented Into the program for the study of numerical characteristics of the 
various constitutive models. Since VPMODEL was written on a modular basis, 
addition of any new viscoplastic model or modification of an existing model 
can be made without major coding effort. For this reason, the program serves 
as an efficient tool for development and refinement of viscoplastic models. 

It can also be used for understanding the physical nature of a given consti- 
tutive model as well as evaluating and comparing different viscoplastic models 
and numerical Integration schemes for use In structural analysis for aerospace 
and other Industry applications. 
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APPENDIX A - USER'S INSTRUCTION 


Summary of Analysis Features 

The current version of VPMODEL consists of the following analysis features: 

(a) Miller, Walker, KRS, and Robinson viscoplastic models. 

(b) Uniaxial responses for monotonic or cyclic loading, creep stress 
relaxation or any combination of these phenomena. 

(c) Isothermal or nonisothermal (e.g., thermal -mechanical cycling) 
condition. 

(d) Integration schemes including forward Euler method, explicit 
trapezoidal method, and implicit trapezoidal method with 
Newton-Raphson iteration. 

(e) Conversational data input with free-field format. 

(f) Interactive computations 

(g) Plotting capability on graphics terminal. 


Nomenclature 

In the conversational questions for data input to VPMODEL, several 
terminologies are referenced and they are described as follows: 


Execution - Each time when the user loads the VPMODEL program into 

computer to perform one or several analyses, it is called 
an 'execution'. This is also equivalent to a terminal 
session. 


Case 


Segment 
Hold Time 


Time Limit 


During each execution (or terminal session), the user may 
wish to perform one or several cases of analysis; for 
example. Case 1: Cyclic stress-strain response using 

Walker's model. Case 2: Cyclic stress-strain response 

using Miller's model, etc. 

Each case may consist of one or more load segments. 

The time period during which creep under constant stress or 
stress relaxation under constant strain is being exercised. 

The maximum time at which the analysis of a load segment 
terminates . 


Stress Limit- The maximum or minimum stress value (including the 

algebraic sign) at which the analysis of a load segment 
terminates . 

Strain Limit- The maximum or strain value (including the algebraic sign) 
at which the analysis of a load segment terminates. 
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Limitations 


VPMODEL has several limitations that the user must be aware of: 

(a) Material constants are available only for Hastelloy-X - Walker's 
model. Miller's model and KSR's model and for 2-1/4 Cr-1 Mo - Steel 
Robinson's model. 

(b) The maximum number of data points Is limited to '4000' for the 

variables: a, E, AE, a, k, T, and t. This limitation can be 

removed by changing the array sizes In COMMON/STAVAR/ In the computer 
program. 

(c) Calcomp plotting routines or Interface routines for IBM TSS/370 
graphics system are required. 


Instructions 

Since the data Input procedure for VPMODEL was written In a conversational 
form, no additional Input Instructions are necessary In order for the user to 
run the program. For each terminal session, the user needs only to provide the 
analysis data In answering those questions prompted by the program and which 
appear as a menu on the terminal screen. The following steps are followed In 
order to run the program on Lewis' IBM 370-3033. 

Step 1 Logon IBM 370-3033 

Step 2 Load the program Into the system by entering "vpmodel', a procdef 
which compiles the program and prepares It for execution. 

Step 3 Provide Input data by answering the questions prompted and 
execute the program for the case being analyzed. 

Step 4 Provide Input data for obtaining results (graphics, listings, 
hard copy, etc.) 

Step 5 Terminate the analysis (or execution) by entering 'O' (for stop) 
from procedure menu. 

Step 6 Logoff 
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TABLE I. - DIFFERENT TEST CONDITIONS FOR 
MODEL COMPARISON 









1 


2 (Load Segment) 




Figure 1. - Different loading histories in creep tests. 




A New Case 



Figure 2. - A flow chart for VPMODEL 
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Figure 2. - Concluded 






STRESS (PSl)-10 



STRAIN -lO-» 


Figure 3. - Stress-strain response of Hastelloy-X, 
Walkers's model: e ■ 3o87 x 10"3/sec; T ■ 

1600° F (8710 C). 
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Figure 4. - Back stress vs inelastic strain for 
Hastelloy-X, Walker’s model: £ • 3.87 x 10“^/ 
see; T = 1600O F (871° C). 
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Figure 5. - Back stress vs strain for Haste! loy-X, 
Walker's model: i - 3.87 x 10“%ec; T • I 6 OOP F 
(871° C). 
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Figure 6. - Cyclic stress-strain response of 
Hastelloy-X, Walker’s model: T ■ 160CP F 
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Figure 7. - Strain history for the stress relaxation test 
of Hastelloy-X; T * 1600^ F (871® C); i ■ 3.87 x 10“^ 
sec. 
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Figure 9. - Stress-strain response 
of 2-1/4 Cr-1 Mo steel, Robinson's 
model e • 2.4 X 10'2/hr, T - 1000“ F 
(538° C). 
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Figure 10. - Back stress vs inelastic strain, Robin- 
son's model: £ ■ 2.4 X 10“ ^/hr, T - 1000^ F 
(538° C). 
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Figure 11. - Cyclic response of 
Robinson's model: T ■ 1000^ F 
(538° C). 
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Figure 12. - Creep response of 2-1/4 Cr-1 Mo steel, Robinson's 
model : e = 0. 6 x 10" ^/hr, T ■ lOOQO F (538® C). 
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Figure 13. - Stress-strain response of 
Hasteiioy-X: T ■ 1400P F (760“ C), 
e ■ 3. 87 X 10'3/sec. 
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Figure 14. - Stress-strain response of 
Hastellov-X: T - 1400° F (760° C), i - 
3. 66 X 10"^/sec. 
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Figure 15. - Stress-strain response of 
Hastelloy-X: T ■ 1400“ F (760P C). i - 
0. 125 X 10'5/sec. 
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Figure 16. - 
Hastelloy- 
3.87x10“ 


Stress-strain response of 
X: T - 1800O F (9820 c), £ - 
^/sec. 
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Figure 17. - Stress-strain response of 
Hastelloy-X: T - 180QO F (982° C), £ ■ 
3. 66 X 10-4/sec. 





Figure 18. - Comparison of step sizes using 
forward Euler method: Walker's model 
T • 1400° F (760“ C), £ - 3.87 x lO'^/sec. 
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Figure 19. - Comparison of different integration 
schemes; forward Euler, explicit and implicit 
trapezoidal methods. 
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Figure 21. - Strain history for dosed symmetric ioop. 
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Figure 22. - Strain-temperature history for closed symmetric loop. 



0*9 



»oi>* (isd) ssayis 


I 


Figure 23. - Predicted response using Walker's model In comparison 
with experimental result for closed symmetric loop. 




Figure 24. - Predicted responses using Miller's model in comparison 
with experimental result for closed symmetric loop. 




Figure 25. - Predicted responses using KSR's model in comparison 
with experimental result for closed symmetric loop. 
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Figure 28. - Strain-temperature history for open symmetric history. 
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Figure 29, - Predicted response using Walker's model in com- 
parison with experimental result for open symmetric loop. 
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Figure 30. - Predicted response using Miller's model in com- 
parison with experimental result for open symmetric loop. 
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Figure 31. - Predicted response using KSR's model in compar- 
ison with experimental result for open symmetric loop. 




Figure 32. Temperature history for open non-symmetric loop. 
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Figure 33. - Strain history for open non -symmetric loop. 




Figure 34, - Strain-temperature history for open non- 
symmetric loop. 
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Figure 35. - Predicted response using Walker's model in com- 
parison with experimental result for open non-symmetric loop. 
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Figure 36. - Predicted response using Miller's model In com- 
parison with experimental result for open non-symmetric loop. 




Figure 37. - Predicted response using KSR's model in compari- 
son with experimental result for open non -symmetric loop. 
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